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ABSTRACT 

We investigate the dynamics and morphology of jets propagating into the 
interstellar medium using 2D relativistic hydrodynamic simulations. The calcu- 
lations are performed assuming axisymmetric geometry and trace a long propa- 
gation of jets. The jets are assumed to be 'light' with the density ratio between 
the beam to the ambient gas much less than unity. We examine the mechanism 
for the appearance of vortices at the head of jets in the hot spot. Such vortices 
are known as a trigger of a deceleration phase which appears after a short phase 
in which one dimensional analysis is dominant. We find that an oblique shock at 
the boundary rim near the end of the beam strongly affects the flow structure in 
and around the hot spot. Weakly shocked gas passed through this oblique shock 
and becomes a trigger for the generation of vortices. We also find the parameter 
dependence of these effects for the propagation and dynamics of the jets. The 
jet with slower propagation velocity is weakly pinched, has large vortices and 
shows very complex structure at the head of the jets and extended synchrotron 
emissivity. 

Subject headings: galaxies: jets — hydrodynamics - shock waves - methods: 
numerical — relativity 



1. INTRODUCTION 

It is widely known that there are three classes of highly collimated and supersonic jets 
from dense central objects with accretion disks. Depending on the central object the following 
can be distinguished : protostars, binary stars, and Active Galactic Nuclei (AGN). AGN jets 
are the largest scale phenomena and the velocity of the jet beam is highly relativistic, at 
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least close to the central objects (see for example, Pearson et al. (1981); Biretta, Sparks, 
& Macchetto (1999)). The jet, which originates near an accretion disk that surrounds an 
AGN, can propagate over a long distance up to a few Mpc while remaining well collimated. 
There are two shocks at the end of the jet. One is a bow shock (or a forward shock) which 
accelerates the ambient gas. The other is a terminal Mach shock (or a reverse shock) at 
which the beam ends. At the terminal Mach shock, non-thermal particles are accelerated 
and emit photons due to synchrotron radiation and inverse Compton scattering. The gas 
which crosses the terminal Mach shock into a hot spot is hottest and pressurized and then 
expands laterally and envelops the beam with the shocked ambient gas, creating a so-called 
cocoon structure. At the contact discontinuity between the ambient gas and the jet in the 
cocoon Kelvin- Helmholtz instabihties develop. 

Cygnus A is a suitable object to see these features because it is one of the closest 
radio galaxies and the beam propagates to perpendicular direction to line of sight. Prom 
observations, its size is about 120kpc, the beam velocity is 0.4 ~ Ic and the hot spot's ram 
pressure advance speed is 0.03c (Carilli et al. 1998), where c is the speed of light. Por more 
details of AGN jets, see for example Carilh & Barthel (1996); Perrari (1998); Livio (1999) 
and references therein. 

Other active regions in jets are knots and secondary hot spots. Some knots with non- 
thermal emission can be seen sporadically along the straight beam flow. The emissivity 
between the knots varies even if they belong to the same jet. Recent high resolution ob- 
servations of AGN jets show fine structure of knots in the beam flow for up to several tens 
of kilo parsecs in M87 (Biretta, Sparks, & Macchetto 1999; Marshall et al. 2002; Wilson & 
Yang 2002), 3C273 (Bahcall et al. 1995; Marshall et al. 2001; Sambruna et al. 2002; Jester 
et al. 2002), the jet in Centaurus A (Kraft et al. 2002), 3C 303 (Kataoka et al. 2003), and 
others (Sambruna et al. 2001). Knots on sub-parsec or parsec scales are thought to be due 
to intermittent flows from the central source and knots in 'blazars' are due to exceptionally 
strong shocks which are caused by the collisions of internal shocks (Spada, Ghisellini, Laz- 
zati, & Celotti 2001; Bicknell & Wagner 2002). Knots on kilo parsec or larger scales are not 
well understood. The high energy particles that are accelerated on parsec scales may retain 
their energy for a long time. In some jets, one observes secondary hot spots adjacent to the 
primary hot spot at the head of jets. The emissivity of secondary hot spot can be as high as 
that of primary hot spot. It is thought that the gas in the secondary hot spot is separating 
from the primary hot spot. Por example, at both eastern and western side of Cygnus A jet, a 
bright secondary hot spot is observed (Wilson & Yang 2000). The reason for these multiple 
hot spots is not well understood. 

Recent observations of AGN jets provide information not only about the features of 
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large scale jets but also that of very small scale jets such as in 'Compact Symmetric Objects' 
(CSOs) which are believed in the first stage of AGN jets from its morphology. CSOs show 
two sided emissivity within the central few kpc of its parent galaxy. The physical properties 
of two hot spots and lobes at both sides have been measured (Owsianik & Conway 1998; 
Polatidis & Conway 2003; Giroletti et al. 2003). CSOs are only a few thousand years old. 
and the propagation velocity of the hot spot is a few tens of percent of the speed of light. 
CSOs may tell us about the physical conditions during the earliest phases of jets. 

To understand the morphology and the dynamics of the jet, analytical studies and 
numerical simulations have been performed for the past thirty years. Blandford & Rccs 
(1974) and Scheuer (1974) discussed the structure of the jet with the relativistic beam 
model theoretically. Some of the early numerical work was done by Norman et al. (1982). 
Although these were non- relativistic simulations, they showed the main structure of jets, that 
is two shocks at the head of the jet and the cocoon. The difficulty of numerical relativistic 
hydrodynamics has delayed the investigation of the relativistic effects on the morphology and 
the dynamics of jets. Only since the past ten years stable codes with or without external 
magnetic fields have been developed for the ultra-relativistic regime (Eulderink & Mellema 
1995; Duncan & Hughes 1994; Font et al. 1994; Koide, Nishikawa, & Mutel 1996; Komissarov 
1999; Aloy et al. 1999a; Hughes, Miller, & Duncan 2002). which allowed the investigations 
of the formation, coUimation, and propagation of jets (Rosen et al. 1999; Koide et al. 2002) 

Marti et al. (1997) and Aloy et al. (1999b) performed long term simulations to study 
the morphology of the jets with 2D and 3D relativistic hydrodynamic codes. Propagation 
of the relativistic jets in the external magnetic field were computed by Komissarov (1999). 
These studies verified that the morphology of jets depends on a number of dimensionless 
parameters: the density ratio of the jet beam to the ambient gas {rj = ph/pa), Lorentz 
factor of the beam, and the Mach number of the beam (Mj, — Vb/cg^b)- Here the subscripts 
'b' and 'a' stand for beam and ambient gas, respectively. In most of these simulations, a 
pressure matched jet is assumed, with the pressure ratio K = Pb/pa = 1- Marti, Miiller, 
& Ibariez (1998) performed long term simulations of jets and found a deceleration phase of 
the propagating jets. Recently Scheck et al. (2002) have pointed out that the propagation 
velocity, Vj^, which is derived from one dimensional momentum balance at the rest frame 
of the contact discontinuity is very useful for understanding the propagation of the jet into 
the ambient gas. They assumed a constant velocity Vj^ — 0.2c and the kinetic luminosity 
of the beam is fixed to be Lkin — lO^^erg s~^ in their models. Those are the first numerical 
simulations using relativistic hydrodynamic code which assume the propagation efficiency 
Vj'^/vb is less than 0.5 (cf. most numerical simulations done so far asuume the efficiency is 
0.5 ~ 1). Thus more wide range of parameter space of vj^ is necessary if we agree with 
recent observations of CSOs are young AGN jets. Scheck et al. (2002) also investigated 
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the dependence of the composition of jets using a generic equation of state. They varied 
the density ratio (rj), pressure ratio (K), Mach number (Mb) and the fraction of protons, 
electrons, and positrons in their three models. The main difference between the various 
models is seen in the radiative properties of the jets. The deceleration phase is also confirmed 
and shown to begin with the generation of a large vortex at the head of jet. It seems 
indeed consistent with observations that jets have a deceleration phase because CSOs have 
a propagation speed of a few tens percent of the speed of light, while in large scale jets such 
as Cygnus A the velocity is only a few percent of speed of light. However the physics of the 
generation of vortices is not well understood. 

A first attention to compare the vortex structures in numerical simulations of the jets 
with observational data have been made by Saxton et al. (2002); Saxton, Bicknell, & Suther- 
land (2003). who discussed the shock structures at the head of the jet of Pictor A and rings 
in the cocoon of Hercules A. We note that radio jets with such features are rare (Gizani & 
Leahy 2003). 

In this paper we will examine the mechanism that causes such vortices in some detail. 
We present three simulations with different vj^ as injected beam conditions. We pay special 
attention to the generation of vortices in the hot spot and the separation of vortices from 
hot spots which strongly affects the dynamics of jets. We show the mechanism of generation 
of vortices at the head of jet. We hope that this study will help us to better understand the 
propagation of jets into the ambient gas. This paper is organized as follows. In Sect. 2, the 
basic equations arc presented. In Sect. 3, we show our three models for the numerical nu- 
merical simulations. The results and discussion are shown in sect. 4. The surface brightness 
of synchrotron emissivity and the resolution study are also shown. The conclusion is given 
in Sect. 5. 



2. BASIC EQUATIONS 

2.1. RELATIVISTIC HYDRODYNAMIC EQUATIONS 

Relativistic flows, ~ 1, affect the morphology and dynamics of jets. The relativistic 
hydrodynamic equations (Landau & Lifshitz 1987) then have to be solved: 
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d(phW^-p) ^ I dr(phW^Vr) ^ djphWv,) _^ 
dt r dr dz ^ 

where, p,p,Vi, W, and h are rest mass density, pressure, three velocity component in i direc- 
tion, Lorentz factor (= (1 — v"^)^^^"^), and specific enthalpy (= 1 + e + p/p), respectively. In 
this section we equate the velocity of light to unity. We assume axisymmetry. The set of 
equations (1-4) is not closed until the equation of state is given. In this paper we assume an 
ideal gas; 

P = (7 - (5) 
where 7 and e are a constant specific heat ratio and specific internal energy, respectively. 



2.2. PROPAGATION VELOCITY 

A one dimensional analysis is useful for estimating the propagation velocity of the jet. 
Here the propagation velocity Vj is defined to be the velocity of the contact discontinuity 
at the end of the jet in the rest frame of the ambient gas. For very 'fight' jets {rj <^ 1), in 
general, the propagation velocity is less than that of the beam because the conversion of the 
large fraction of the beam kinetic energy to the thermal energy occurs at the terminal Mach 
shock. 



2.2.1. NON RELATIVISTIC CASE 

If the velocity and temperature of the gas are in the non-relativistic regime, it is sufficient 
to consider the Euler equation. Assuming momentum balance in the rest frame of contact 
discontinuity, 

Sb {pbiVb - v/^'f+Pb) = Sa {pa{-v/''f+Pa) , (6) 

where Sb and Sa are the cross section of the beam and the hot spot, and t;/^ is the propa- 
gation velocity derived from one dimensional analysis. Solving Eq.(6) for w/^, we obtain 

Am. jAhfvl {Ar, l){Ajrvl + (AK mjj) 

= A^l ' 

where Cg is sound speed and A is cross section ratio {A = Sb/Sa)- If we can assume A—1 
and neglect the term that includes K{= Pb/Pa) in Eq.(6), then this equation reduces to 
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This is the same as the equation derived by Norman, Winkler, & Smarr (1983) for the 
pressure matched jet {K =1). If the density ratio r] is much less than unity, in the case of 
'light' jet, the propagation efficiency defined as Vj/vb is much less than unity. 



2.2.2. RELATIVISTIC CASE 

For the relativistic case, assuming again ID momentum balance between the beam and 
the ambient gas in the rest fame of the contact discontinuity at the head of the jet, 

Sb {pbhbWfW^ivb - vf f + pb) = Sa {pahaWfi-vf f + p„) , (9) 

where Wj = (1 — vj^^)"^''^. This equation leads to 
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In general, the sound speed of interstellar matter is much smaller than the velocity of the 
beam (cg^a <^ "^b) so we can neglect the term including K. As a result, Eq.(lO) becomes 



vf = -^^^Vb. (12) 



This is equal to the equation derived by Marti et al. (1997) for the pressure matched jet 
[K = 1) if we assume A = 1. Scheck et al. (2002) showed that the time evolution of the 
actual propagation velocity is almost the same in their models for fixed Vj^ and A = 1, 
in spite of the large variation of the specific internal energy of the beam during the first 
phase. There are several ways to explain the deceleration of jets between sub parsec and 
mega parsec scales. If the density ratio rj varies, rj should become smaller than that in the 
earher phase. The Lorentz factor should also become smaller. The velocity at the head of 
the beam should be relativistic to reproduce strong emissivity at the hot spot. The effect 
of the Lorentz factor seems small. Another possibility is variation of the density profile of 
the ambient gas. For example, it is natural to assume that the ambient gas decreases in 
density as function of distance from the AGN. In the models of Scheck et al. (2002), the 
multidimensional effect, namely decreasing cross section's ratio, becomes very important for 
the dynamics and the jet decelerates gradually. Scheck et al. (2002) compared the time 
evolution of the dynamics with the extended Begelman-Cioffi model. The original model 
derived by Begelman & Cioffi (1989) assumes a constant propagation velocity and Scheck et 
al. (2002) extended the assumption to a power law. 
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3. NUMERICAL CONDITIONS 

We solve the relativistic hydrodynamic equations Eq.(l-5) with a 2D relativistic hy- 
drodynamic code recently developed by one of the authors (A.M.). The detailed numerical 
method and test calculations of the code are given in the appendix of this paper. 

We use a 2D cylindrical computational region r x z. The grid size in both r and 
z directions is uniform, namely Ar = A2; = const. We assume that the ambient gas is 
homogeneous initially. Calculations including clouds in the ambient gas have been made 
by de Gouveia Dal Pino (1999); Zhang, Koide, & Sakai (1999); Hughes, Miller, & Duncan 
(2002). A relativistic beam flow (wf,) which is parallel to z axis is continuously injected from 
one side of the computational region {z — Q). The inner 10 grid points from the symmetry 
axis are used for this. The radius of the injected beam = lOAr is used as a scaling unit 
in this study. The computational region covers 50R]^(r) x 180R]^(2;); corresponding to a grid 
500 X 1800 grid. The radius of the beam near the central engine is unknown, and we assume 
a plausible value of 1R|^ = O.Skpc, corresponding to 25kpc x 90kpc for the computational 
region. A free outflow condition is employed at the outer boundaries for r = 50R|^ and 
z — ISOR]^. A reflection boundary is imposed on the symmetry axis as well as at 2; = with 
r > R^. The boundary condition at 2; = is crucial for dynamics and the outer shape of 
the jets. The free boundary condition permits the gas to escape at the back side (Saxton et 
al. 2002). According to observations jets have counter jets which propagate in the opposite 
direction. We assume the reflective boundary condition so that our calculations begin near 
central engine. 

We examine three injection beam conditions (Table 1) with different v/^. The de- 
scription of Vj^^ is given in Eq.(12), here A — 1 is assumed. The beam velocity is fixed 
to Vb = 0.99c (the corresponding Lorentz factor is Wb = 7.1). The Mach number of the 
beam is also fixed to Mb = 6.0. We assume that the beams don't have highly relativistic 
temperatures and that remains constant Vb- As a result Vj^^ depends only on the density 
ratio 1]. Scheck et al. (2002) adopted Vj^^ = 0.2c as the injected beam conditions for their 
calculations. This is in the range of acceptable propagation velocity based on observation of 
CSOs. We explore the region of the propagation velocity vj^^ from 0.2c to 0.4c and inves- 
tigate the dependence of the difference of v/^ caused by different 77. With increasing Vj^^, 
rj varies about a factor of ten; from 1.28 x 10"^ to 9.15 x 10~^. We label theses cases JB02, 
JB03, and JB04, and the corresponding Vj^^ are 0.2c, 0.3c, and 0.4c. Although there is only 
about one order difference in the density ratio for all beams, doubling value of vj^^ has a 
big effect on the dynamics of jets and their long term propagation. Because the spatial scale 
is fixed, we adopt different expiration times for each models, see Table 1 again. 

The initial rest mass density of ambient gas Pa is unity for all models. The pressure 
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ratio(ir) is chosen to be from X = 10 to 100 so that the ambient gas is approximately similar 
in all cases. As a result, the temperature of the ambient gas for all injected beams is a few 
eV, if we assume the ambient matter is pure hydrogen gas. 

Because we choose similar values of Vj^ to those of Scheck et al. (2002), other parameters 
(77, K, Mf,) are also very similar. 

4. RESULTS AND DISCUSSION 

4.1. MORPHOLOGY AND DYNAMICS 

Prom our simulations, we find that the overall morphology and dynamics of the jets 
are similar to that discussed in previous work. Figures 1, 2 and 3 show snapshots of rest 
mass density, pressure, and Lorentz factor for three models near the end of the simulations, 
t = 1770 [%/(c)] (model JB02), t = 1140 [%/(c)] (model JB03) and t = 570[%/(c)](model 
JB04). Since these jets have different beam conditions, the end times are also different as 
discussed in the previous section. However this is not only due to different Vj^^ values, 
but also because of differences in deceleration phase of the jets which are caused by the 
generation and separation of vortices at the head of the jet, as we will discuss later. At first 
glance, the outer shapes of the jets are quite different. Model JB02 shows a conical outer 
shape which is very similar to the results found by Scheck et al. (2002). This is not seen in 
the other two models (JB04 and JB03). All beams remain coUimatcd from z = 0, where the 
beam is injected into the computational region, to the head of the jets. It is also important 
to note that the beam radius dose not increase monotonically from the source to the head 
of jets. At the head of the jet the radius of the beam is 3 ~ 5Rb. The opening angle is 
very small 1° ~ 2° for all models. High Lorentz factors exist only along the z axis in the 
beam flow, but they do vary. That means that a large part of the injected kinetic energy is 
transported from the central engine to the head of the jet along the beam flow. At the end 
of the beam a strong 'terminal Mach shock' can be seen. One of the most active points is 
called a 'hot spot' into which shocked beam gas enters through the terminal Mach shock at 
the head of the jet. The pressure in the hot spot is very high due to the energy dissipation 
at the terminal Mach shock and is matched by ambient gas compressed at the bow shock. 
Moderate Lorentz factors {W ~ 2) exist in the cocoon. These are due to back flows that 
begin at a hot spot and flow back in the center of the cocoon parallel to the beam flow. This 
mildly relativistic back flow seems a little strange thing because the head is proceeding very 
slow (< 0.2) and the expanding velocity from the hot spot is up to 0.5c in the comoving 
frame (maximum sound speed). It should be noted that this mildly relativistic backflow 
is longer for slower jets. Thus some acceleration mechanism or other effects are necessary. 
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We discuss this in the next subsection. This back flow creates a shear flow and the contact 
surface becomes unstable due to Kelvin-Helmholtz instabilities. The surface between the 
back flow and shocked ambient gas's flow also becomes unstable and causes the appearance 
of vortices in the larger cocoon. 

To understand the difference of these models after long term propagation we must study 

the dynamics of the jets. Figure 4 shows the time evolution of the positions of the contact 
discontinuity the bow shock (forward shock), and the terminal Mach shock (reverse shock) 
at the head of the jet for each model. Three straight lines which correspond to the lines 
of one dimensional analysis are also shown for the comparison. It is difficult to define the 
positions of these surfaces because of the complex structure at the head of the jets. In this 
paper we show the positions at r = 0, namely along the z axis. It is hard to identify the 
contact discontinuity between shocked ambient gas and shocked beam gas especially in later 
phase because of mixing of the shocked beam with the ambient gas due to the generation of 
vortices. The contact discontinuity is defined as the boundary of shocked ambient gas that 
has half the maximum density. We plot the points where the Lorentz factor becomes 2 at the 
head of the jet. The two phases indicated by Marti, Miiller, & Ibaiiez (1998) and Scheck et 
al. (2002) are also seen in our results. During the first phase, all of the slopes are constant. 
Model JB02 has a propagation velocity ~ 0.2c and the one dimensional analysis discussed 
in section 2.2 is in good agreement with our result. For JB03 and JB04 the propagation 
velocity is a little faster than that of in our one dimensional estimates. During first phase 
the surfaces are very close each other. On the other hand, in the second phase, the surfaces 
separate and approach each other repeatedly. The propagation velocity is no longer constant 
and the jet is decelerating gradually. Some vortices grow at the head of jet and separate 
towards the back side. This affects the dynamics of the jets. 

Figure 5 shows snap shots of the velocity {^Jv^ + vf) during the earlier phase of the 
each models. Although its path is not as straight as that of the beam, the back flow does 
remain parallel to the beam flow. The back flow velocities are very similar for all models 
(0.3 ~ 0.4c). Because the propagation velocity is less than that of back flow for the model 
of JB02, most of it reaches the boundary, after which some of the gas expands in lateral 
direction, and other gas form a new fiow which is between the beam and the back fiow. This 
new flow also can be seen in the case of JB03 near the boundary, but it does not reach the 
head of jet hke JB02 because the propagation velocity is faster than that of JB02. The beam 
flow and back flow he side by side in the model JB03 and JB04 near the head of the jets. 

In later phases this 'third' flow, which is not beam flow but has the same propagating 
direction, exists in all models because of vortices generated in the cocoon due to Kelvin- 
Helmholtz instabilities. At this stage the back flow becomes very complex. How much gas 
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reaches the boundary at r = affects the expansion into lateral directions for each model. 
Slower propagation velocities result in cone like outer shapes of the jets. The shocked ambient 
gas forms a shell-like structure like a football, resembling the observed X-ray emission Cygnus 
A by Chandra (Wilson & Yang 2000). It is consistent with our results because the present 
propagation velocity of Cygnus A is very small (~ 0.03c). A long time may have passed after 
the deceleration phase began in the Cygnus A jet. In contrast, in models JB03 and JB04, 
which have propagation velocities faster than that of the back flow in the earlier phase do 
not have a conical shape. These jets are decelerating. The effect of a decreasing propagation 
velocity can be seen near the head of the jet in JB03, showing a conical shape around the 
head of the jet. 



4.2. VORTEX FORMATION IN HOT SPOTS 

We discussed the generation and separation of large vortices at the head of the jet which 
occurs repeatedly during the second phase. Such processes strongly affect the dynamics of 
jets and was also found by Scheck et al. (2002). It is also consistent with observations between 
the high speed CSO sources and slower large scale jets. The reason for this deceleration seems 
to be the formation of vortices. Where and how are these vortices created ? The possibilities 
of hydrodynamic instabilities at the head of jets have been discussed before. Norman et 
al. (1982) discussed that there are Rayleigh- Taylor instabilities at the contact discontinuity. 
Recently, Krause (2003) mentioned the Kelvin-Helmholtz instability between the flow from 
hot spot and shocked ambient gas. 

To understand the generation of these vortices we need to focus on the flow structure 
in and around the hot spot. The vortices appear in the hot spot and they originate by the 
flow through an oblique shock at the end of the beam. 

It is known that oblique shocks appear within the beam, even before the end, when 
the beam expands in lateral directions. This re-establishes pressure equilibrium between the 
beam and its surroundings and keeps it confined (Fig 6a). When the beam is pinched for 
some reason, the gas tends to expand due to increased pressure caused by the compression 
(Fig 6b). If the pressure outside is high enough to confine the beam, an oblique shock 
appears to prevent it from expansion. (Fig 6c). The radius of the beam after reconfinement 
depends on the pressure outside of the beam and how the beam is pinched. 

Figure 7 shows the rest mass density, pressure and Lorentz factor profile along the z 
axis of model JB02 at t = 300, 600, 900, 1200, and ISOOR^/c. The beam is confined by the 
pressure of the cocoon and shocked ambient gas. However the surface of the beam is not 
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stable. In our calculations the first oblique shock appears where the beam is injected. Such 
shocks appear irregularly in the beam. As these shocks have different speeds, a slower shock 
is caught up by a faster shock. This might explain the high cmissivity knots in the jets. 
When an oblique shock appears at the end of the beam, most of the gas passes a terminal 
Mach shock with large dissipation. However, some of the gas passes through the oblique 
shocks on the side (Fig. 8 (a)-(b)). If the angle is small, the oblique shock is very weak. 
The loss of kinetic energy should be small but the pressure becomes as high as that in the 
cocoon. Then a mildly relativistic back flow begins. The effect of these obhque shocks can 
be seen in slower jets. This is why the mildly relativistic flow is longer for slower jets. Such 
a fast velocity flow through the oblique shock can propagate further in lateral direction than 
slower gas from the hot spot before it propagates backwards. The flow path becomes then a 
circular, arc like vortex. This vortex can sometimes, which triggers instabilities at the beam 
surface and internal oblique shocks, reach the beam flow. The surface of the beam becomes 
more unstable and has oblique shocks there. The oblique shocked flow has an important 
effect on the gas in the hot spot since it blocks further outflow. In the hot spot the velocity 
is not constant. The nearer the gas is to a corner that is a crossing point between a contact 
discontinuity and z axis, the slower the velocity. Although the outflow from the hot spot is 
blocked, the beam flow continues to enter the hot spot. This then drives a clockwise rotation 
vortex in the hot spot (Fig. 8 (c)-(e)). This vortex is very important for the dynamics because 
the increasing radius of vortices causes an increase of cross section of the beam. As a result, 
the beam decelerates (see Eq.(12)). The vortex can grow till the radius becomes about twice 
large as that of the beam. The slower a jet is, the larger vortex it can support because the 
shocked ambient gas which is just going to backward has smaller momentum to this vortex. 
When a vortex grows, the jet decelerates. On the other hand when a vortex separates (Fig. 
8 (f)-(g)), the jet accelerates shghtly. But the next vortex grows soon. In Fig. 4 the size of 
hot spot, namely, the distance between a contact discontinuity and a terminal Mach shock, 
is oscillating with increasing amplitude. This occurs in all of models but is strongest in the 
slower jets. The increase of the amplitude of the oscillation increases the growth time of a 
vortex. A long growth time causes a larger radius of the vortex. The dynamics is a repetition 
of these processes during the second phase. 

Most of gas in this vortex has passed trough the 'terminal Mach shock'. The vortex 
also contains non thermal particles accelerated by this shock. A separating vortex can be 
observed as an active region around the hot spot. The best candidate for this region is 
the secondary hot spot discussed in the introduction. The high energy particles loose their 
energy quickly, so that a separating vortex is a little less bright than that of primary hot 
spot. 

When a jet accelerates, the bow shock has a nose cone like shape. On the other hand. 
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the bow shock which is decelerating has a flat shape. This may be observed in the eastern 
lobe of Cygnus A, where the hot spot seems to be slightly ahead of the lobe. This may be 
because the vortex may just have separated and the jet is in a brief accelerating phase now. 

Once a jet exhibits such oblique shocks and vortices, the flow structure around the hot 
spot becomes more complex. The effect of this turbulence is strong in the case of slowly 
propagating jets because the confinement effect of shocked ambient gas is weak. The beam is 
pinched when a circular arc like flow from the oblique shock reaches the beam or a separated 
vortex moves backwards. The end of the beam tends to expand and has more oblique shocks 
there. The beam is sometimes constricted in the middle by such effects and a strong shock 
appears. This shock may be related to knots observed in the beam. 

During the later phase of the model JB02, the size of hot spot becomes very large and 
clumpy. As a result, the gas which passed the terminal Mach shock follows the edge of this 
clumpy gas. Some of it is mixed in at locations where the gas reaches the shocked ambient 
gas and back flow begins. This then also separates from the head of the jet and 'normal' hot 
spot appears again. The large clumpy gas does not display a bright 'hot spot'. It corresponds 
to a 'temporary absent' hot spot discussed by Saxton, Bicknell, & Sutherland (2003) and 
can be a good candidate to the absent hot spot in Hercules A. 

4.3. SYNCHROTRON EMISSIVITY 

Recent some numerical simulations of jets show the surface brightness of the synchrotron 
emissivity (Scheck et al. 2002; Saxton et al. 2002; Saxton, Bicknell, & Sutherland 2003; Aloy 
et al. 2003). We follow the analysis by Saxton et al. (2002); Saxton, Bicknell, & Sutherland 
(2003) and compare our results with theirs and observations. 

We use an approximation which assumes the synchrotron emissivity is proportional 
to where p is pressure, B is magnetic flled, and a is the spectral index. We use 

a — 0.6 that is typical value. The magnetic pressure ~ is assumed in equilibrium with 
thermal pressure p. The surface brightness is derived from revolved 2D numerical results. 
The emissivity is /p^'^, where / is the fraction of the gas which originats the beam. The 
advection equation of / is transformed to conservative form using mass conservation equation 
and solved with hydrodynamic equations. 

Figure 9 and 10 show the emissivity of the models JB02 and JB04 in log scale with four 
order magnitude from maximum intensity using gray scaled color bar (white is the brightest 
emissivity) in each panel is shown. Different several angles, ^ = 0, 15, 30, 45, 60, 75, and 90°, 
where 9 is the angle between z axis and the line of sight are assumed. JB02 has very extended 



-13- 



emissivity from the head to the root of the jet. On the contrary, JB04 which indicates less 
complex structure has the emissivity only around the head. Thus, the appearance of vortices 
causes the extension of emissive regions toward backside. The brightest region corresponds 
to the hot spot. The other bright regions which corresponds to the oblique shocks in the 
beam appear irregularly. This may correspond to the observed knots in the beam. The same 
emissivity is usually seen in the observations of the jets. Both of models have a ring like 
emissivity near the head of the jet when the angle ^ is 30 ~ 75°. This feature is also shown 
and discussed the similarity with observation of Hercules A in Saxton et al. (2002); Saxton, 
Bicknell, & Sutherland (2003). Several rings in the lobe has different radius and emissivity. 
The observations of Hercules A by Gizani Sz Leahy (2003) resemble our results of emissivity 
in JB02. These results indicate that the lobe of Hercules A has a complex structure with 
vortices. 



4.4. RESOLUTION DEPENDENCE 

At last, we discuss the dependence of the resolution for our discussion. The calculations 
with higher resolution allow to have finer structures because how much numerical dissipation 
is included depends on the resolution. Especially for the problem which has the appearance 
of vortices, such an effect is very important and may affect the size of vortices and the 
timing of formations of vortices. We performed calculations with 1.5 times and twice higher 
resolution for our three cases discussed in the previous sections, but the scales are restricted 
to save CPU time. In other words 15 (1.5 times resolution) and 20 (twice resolution) grid 
points are used for innerlet beam. 

Higher resolution calculations also show the same property as discussed in the previous 
subsections, namely, two phases are observed, the jets have the generation and separation 
vortices at the head of the jet and these effects are strongly seen in the slower injected beam 
models. Figure 11 shows the time evolution of the positions of the bow shocks, contact 
surfaces, and terminal Mach shocks with different resolutions for the condition of JB02. 
This figure is as same as Fig. 4, but only JB02 case is shown and the scales are restricted. 
The jets with higher resolution propagate a little slower than that of normal version. This 
is caused by a large separation of a terminal Mach disk from the contact surface in the early 
second phase because the large separation allows the appearance of large hot spot which is 
seen in the later phase of the calculation in JB02 with normal resolution and large cross 
section for the efficiency of the propagation. 
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5. CONCLUSION 

Three numerical simulations of relativistic jets are shown in this paper. We pay special 
attention to the formation of vortices and its parameter dependence. The propagation ve- 
locity derived from one dimensional momentum balance is varied from 0.2c to 0.4c. These 
estimations for the velocity are based on recent observations of CSOs. We identify two phases 
based long term simulations, confirming previous results. The propagation velocity during 
the first phase may be estimated using a one dimensional analysis. The length between the 
bow shock and the contact discontinuity does not change significantly. On the other hand, 
the distance between the terminal Mach shock and the contact discontinuity oscillates and 
its amplitude increases gradually during the second phase. This effect is strongest in slow 
jets. During the second phase, an oblique shock with a small angle at the end of the beam 
appears. This flow dramatically affects the dynamics and morphlogy, and emissivity of the 
jets in the second phase. Weakly shocked gas then foams a back flow without large energy 
dissipation. The velocity of this flow is faster than that of expanding flow from the hot spot. 
The hot spot is blocked by this weakly shocked flow. As a result, a vortex forms in the 
hot spot. When a vortex grows, the increasing cross section decelerates the jet. When a 
vortex separates from the jet, the jet accelerates. After this a new vortex appears and this 
process occurs repeatedly during the second phase. The jet decelerates gradually in time. 
As a future work, the effect of the generation and separation of vortices at the head of the 
jet should be investigated using three dimensional numerical simulations. 

The synchrotron emissivity of the jets are shown The strong emissivity of the syn- 
chrotron radiation at the hot spot appear. The irregular emissivity on the z axis appears 
and had same similarity with observed knots in the beam. In the slowest beam model very 
extended emissivity is shown from the head to the root of the jet and the ring like structure 
of unusual emissivity of Hercules A is reproduced. In the fastest beam model the emissivity 
can be seen only around the head of the jet and along the beam. The separation of vortices 
from the head of the jet strongly affects the extension of the emissivity. 

This work was carried out on NEC SX5, Cybermedia Center and Institute of Laser 
Engineering, Osaka University. We appreciate computational administrators for technical 
supports. 
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In this appendix, we describe a numerical method used in our relativistic hydrodynamic 
code and show some results of typical test problems for relativistic hydrodynamic codes, 
namely, the shock tube problem and strong reflection shock problem. 

During the past ten years, numerical methods to solve the relativistic or magneto- 
relativistic hydrodynamic equations have advanced significantly (see a review paper by 
Ibaiiez & Marti (1999) and references therein). Especially, the method using approximate 
Riemann solvers provides good accuracy even if the flow includes strong shocks and high 
Lorentz factors. Recently we have developed a two-dimensional special relativistic code using 
approximate Riemann solvers which are derived by spectral composition of Jacobian matrices 
of special relativistic hydrodynamic equations. Plane, cylindrical(r — ^), and spherical (r — ^) 
geometry are assumed. 

Recently other new methods also have been employed to relativistic hydrodynamic or 
relativistic magnet ohydro dynamic codes and proposed (Sokolov (2001); Del Zanna & Buc- 
ciantini (2002); Anninos & Fragile (2003); Del Zanna, Bucciantini, & Londrillo (2003)). 



Relativistic hydrodynamics equations are written in conservative form in each coordinate 
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where, w, / and gf, and Sc and Sg are conservative vector, flux vectors, and source vectors, 
respectively. They are defined as follows, 

u = {pW, phW\\ phW\^, phW^ -p- pWf, (A4) 

f{u) = {pWv\ phWW^ + p, phW^vW, phW\^ - pWvY, (A5) 

g{u) = {pWv^, phW\h^, phWW^ + p, phW^v"^ - pWvY, (A6) 

Se= (0,^,0,0)^, (A7) 

/ 2p + pWVv2 phW^v^v^ -cot Op \^ 



Then, we disperse these equations for the numerical calculations. 



— (n+l/2,i/i+l/2,j -ri-l/2,jfi-l/2,j) -T— 
' i,j ^ ' i 

At 



(A9) 



where m"^ stands for the average value of u in the i,j-th grid at and / and g stand 
for numerical flux through the cell surface. Here, we show the cylindrical coordinate case. 

Numerical fluxes / and g at each cell surface are calculated by Marquina's flux formula 
derived from left and right eigenvectors and eigenvalues of Jacobian Matrices of relativistic 
hydrodynamic equations (Donat & Marquina 1996; Donat et al. 1998). To obtain higher 
accuracy in spatial dimensions, we adopt the MUSCL method (van Leer 1977, 1979) for the 
reconstruction of the left and right state at each cell surface. 

{qL)i+i/2 = + ^ ((1 - «)(A_), + (1 + K){K+)i) , (AlO) 

{qR)i+i/2 = Qi+i - ^ ((1 - «^)(A+)«+i + (1 + K){A.)i+i) , (All) 

where q is the physical value for interpolation. In our code the primitive values such as rest 
mass density, pressure, and velocity are used for this reconstruction. The accuracy of the 
code for spatial dimensions is second order when a linear interpolation is adopted {k = —1 
or, 0), and third order when quadratic function is used for the interpolation (k = 1/3). We 
use minmod limitcr to keep Total Variation Diminishing (TVD) condition which prevents 
the development of numerical oscillations. Then A_|_ and A_ are defined; 

(A+)j = minmod(5j+i - g^, h{qi - ?i_i)), (A12) 
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(A_)i = minmod(gi - g^-i, b{qi+i - Qi)), (A13) 
mmmod(a, 6) = sign(a)max(0, min(|a|, sign(a) b)), (^14) 

where 6 is a parameter that satisfies 1 < < (3 — — k). The resuhs become diffusive 

with small b. In this study we use k = — 1 and 6 = 2. The accuracy of the time-step is first 
order. Recovery of primitive values from conservative vector u is done by Newton-Raphson 
method at every time step (Aloy et al. 1999a). 

A.2. TEST CALCULATIONS 

We show two types of test calculations of 1-dimensional relativistic hydrodynamics. The 
calculation is done with a version of 2nd order for space and 1st order for time. 400 grid 
points in the interested direction are used for both cases. As we assume the dynamics is 
only one dimensional, the velocity in the other direction is set to be zero initially. So, the 
velocity of this direction remains zero with time evolution. 

A.2.1. SHOCK TUBE PROBLEM 

A shock tube problem is a kind of initial value problem. Two states are given between 
a discontinuity ai t — initially . This is reasonable because the analytical solution is given 
by Thompson (1986); Marti & Miiller (1994). We use two initial conditions of this problem 
following Donat et al. (1998). 

• Problem 1 : SHI (Plane) 
Left state {x < 0.5); pl 
Right state {x > 0.5); pr 

• Problem 2 : SH2 (Plane) 
Left state {x < 0.5); pl 
Right state {x > 0.5); pr 



\,Pl = 13.3, = 0,7l = 5/3 
l,Pfl = lxlO-6,^;« = 0,7^ = 5/3 



10 



= l,PL = 1000,T;i = 0,7i = 5/3 
-- 1,p^ = 0.01,^;r = 0,7r = 5/3 



Figures 12 and 13 show the results of these problems at t = 0.5 (SHI) and t = 0.35 
(SH2). Analytical solutions are also shown for the comparison. Despite of more than 5 
orders jump in pressure at the initial discontinuity, a rarefaction fan whose front proceeds 
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towards the left with a sound velocity of the left side state of initial condition, a contact 
discontinuity and a shock are given with good accuracy in each case except at the density 
jump behind the shock in SH2. 



A.2.2. REFLECTION SHOCK PROBLEM 

A reflection shock problem is suitable for studying a flow with a strong shock. Initially a 
relativistic homogeneous cold flow (po, W^o, ~ 0) reflects at x = 0; a wall(plane geometry), 
symmetric axis (cylindrical geometry) or symmetric point (spherical geometry). Then, a 
strong shock runs against the flow with a high density and pressure jump. Analytical solu- 
tions of this problem are given by Rankine-Hugoniot relation of relativistic hydrodynamics 
(Johnson & McKee 1971). 

{^^ + ^{W,-l)]po.e = Wo-l,v = Q, iorxKVst (A15) 
\7 - 1 7-1 J 

p = po ( 1 + -^^^ ) , e ~ 0, = ■uo, for X > Vgt, (A16) 



X 



where Vg is shock velocity. 

The geometry is plane for (a = 0), cylindrical (a = 1), and spherical (a = 2). The expression 
of rest mass density jump at the shock front is divided in two parts, namely a maximum 
density compression ratio in non relativistic hydrodynamics ((7 + 1)/ (7 — 1)) and including 
a Lorentz factor. 

• Problem 3 : REP 

po = 1.0, eo = 10-^, vq = -0.999(H/o = 22), 7 = 4/3, Plane 

• Problem 4 : REC 

po = 1.0, eo = 10-^, Vq = -0.999(Wo = 22), 7 = 4/3, Cylindrical 

• Problem 5 : RES 

po = 1.0, eo = 10-^^;o = -0.999(M/o = 22), 7 = 4/3, Spherical 



Figures 14, 15, and 16 show the numerical results for rest mass density, pressure, and 
velocity for each of the problems. Analytical solutions are also given. We show the results at 
t — 1.57 when the shock propagates 0.5 from x — Q. In all cases the shock front is captured 
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with good accuracy. Because of numerical problem, a oscillation appears at the front. At 
the boundary, the error is within 2 % (REP), 6 % (REC), and 5 % (RES) in density The 
geometry error in cylindrical and spherical case is 1 % (REC), and 2 % (RES) in density. 
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expiration time Yly^/{c) 
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Table 1: numerical conditions of models 
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Fig. 1. — Contours of rest mass density (top), pressure(middle) , and Lorentz factor (bottom) 
of model JB02 at the end of the simulation {t = 1770 [R^/c]) 
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Fig. 2. — Contours of rest mass density (top), pressure(middle) , and Lorentz factor(bottom) 
of model JB03 at the end of the simulation {t — 1140[Rb/c]) 
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Fig. 3. — Contours of rest mass density (top), pressure(middle) , and Lorentz factor(bottom) 
of model JB03 at the end of the simulation {t — 570[Rb/c]) 
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Fig. 4. — Time evolution of the position of the bow shock, the contact discontinuity, and 
the terminal Mach shock for models JB02, JB03 and JB04. 
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Fig. 5. — The absolute velocity (r > 0) and log scale rest mass density contour 
(r < 0) of model JB02 (top), JB03 (middle), and JB04 (bottom) in early phase, at 
t = 107.5Rb/c(JB02), t = 50.0%/c(JB03), t = 35.0Rb/c(JB04) 
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Fig. 6. — A schematic of appearance of an oblique shock in the beam (not scaled) 
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Fig. 7. — The log scale rest mass density (solid), log scale pressure (dashed), 
and Lorentz factor (dashed point) profile along the z axis of model JB02 at t= 
300, 600, 900, 1200, andl500Rb/(c). 
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Fig. 8. — A series of profiles of growing and separating of vortices. In each panel the top 
shows the absolute velocity r > and the bottom shows rest mass density r < 0, (a):normal 
profile, (b):an oblique shock appears, (c)-(e):a vortex grows up a vortex, (f)-(g):separates 
from top of the jet 



-31 - 



Fig. 9. — Synchrotron emissivity of JB02 at t = ITTOfR^/c] with several angle 9 = 
0, 15, 30, 45, 60, 75, and 90°, where 9 is the angle between z axis and the line of sight. Coun- 
ters are shown in log scale with four order magnitude from maximum intensity using gray 
scaled color bar (white is the highest emissivity) in each panel. 
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Fig. 10.— Same as Fig. 9 of JB04 at t = STOfRb/c]. 
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Fig. 11. — Time evolution of the surfaces at the head of the jets (same as 4 but only JB02) 
using different resolution for the calculation. A (innerlet beam 10:normal case), B (15), 
C(20). 
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Fig. 12. — Test calculation of shock tube problem (SHI), rest mass density (p), pressure (p), 
velocity (v) are shown at t=0.5. Initial discontinuity is at x = 0.5. Solid lines are analytical 
solutions. 
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Fig. 13.— Same as Fig. 12 of SH2 at t = 0.35. 
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Fig. 14. — Test calculation of reflection shock problem (REP;plane), rest mass density (p), 
pressure (p), velocity {v) are shown at t=1.57. Solid lines are analytical solutions. 
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Fig. 15.— Same as Fig. 14 of REC(cyhndrical) at t = 1.57. 
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Fig. 16. — Same as Fig. 14 of RES (spherical) &X, t — 1.57. 
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